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We discuss tlie effects of quenclied disorder on a model of cliarge density wave (CDW) ordering on 
tlie square lattice. Our model may be applicable to the cuprate superconductors, where a random 
electrostatic potential exists in the Cu02 planes as a result of the presence of charged dopants. We 
argue that the presence of a random potential can affect the unidirectionality of the CDW order, 
characterized by an Ising order parameter. Coupling to a unidirectional CDW, the random potential 
can lead to the formation of domains with 90 degree relative orientation, thus tending to restore 
the rotational symmetry of the underlying lattice. We find that the correlation length of the Ising 
order can be significantly larger than the CDW correlation length. For a checkerboard CDW on the 
other hand, disorder generates spatial anisotropics on short length scales and thus some degree of 
unidirectionality. We quantify these disorder effects and suggest new techniques for analyzing the 
local density of states (LDOS) data measured in scanning tunneling microscopy experiments. 



I. INTRODUCTION 



One of the major stumbling blocks preventing a quan- 
titative confrontation between theory and experiment in 
the cuprate superconductors is the influence of quenched 
disorder on the experimental observations. The dopant 
ions exert a significant electrostatic potential on the 
Cu02 plane, and so unless the ions can be carefully 
arranged in a regular pattern, the mobile charge car- 
riers experience a random potential. Recent STM 
observationiiSi^iii^ clearly display that quenched random- 
ness is crucial in determining the spatial modulations of 
the local density of states. 

There has been much recent interest in determin- 
ing the nature of the spin and charge density wave 
order (CDW) observed in STM, neutron, and X-ray 
scattering in a variety of cuprate compounds at low 
temperaturesiiSiSi^iSiSiLSiS. The quenched disorder acts 
on the CDW order as a "random field", which is al- 
ways a relevant perturbation at low temperatures: true 
long-range order is disrupted at any finite random field 
strengtbiSi. Nevertheless, one might hope that an an- 
alytic treatment may be possible in the limit of weak 
random fields. Many such analysesiSiiiiiSiiii^ have been 
carried out in the literature, describing states with power- 
law correlations and suppressed dislocations (or related 
topological defects) at intermediate length scales. At the 
longest scales, dislocations always proliferate and all cor- 
relations are expected to decay exponentially; no ana- 
lytic treatment is possible in this strong coupling regime. 
As we will discuss below, current experiments on the 
cuprates are in a regime dominated by dislocations, and 
there does not appear to be any significant regime of 
applicability of the defect-free theory. Consequently, we 
are forced to rely on numerical simulations for an un- 
derstanding of experiments. We will present numerical 
results over a representative range of parameters. Our 
aim is to allow insights into the underlying theory by a 
comparison of experimental and numerical results. 



II. MODEL 



A previous work by two of the authorsl^ studied the 
influence of thermal fluctuations on density wave order 
on the square lattice. Here, we will study the influence 
of quenched randomness on the same underlying theory. 
A generic density was defined which could be any ob- 
servable invariant under spin rotations and time reversal 



5p{r) = Re [^^e^'^^-'-] + Re ['^ye'^^'''] , (2.1) 



where K^; = (27r/a)(l/p, 0), Ky = (27r/a)(0, 1/p), and 
'^x,y are complex order parameters which were assumed 
to vary slowly on the scale of a lattice spacing. 

If both amplitudes have nonzero expectation val- 

ues, the charge density is modulated in both x- and y- 
directions and describes a solid on the square lattice. In 
addition, if the wave length of the charge ordering is com- 
mensurate with the underlying crystal, i.e. for integer p, 
the density displays true long range order. For incom- 
mensurate charge order, fluctuations due to flnite tem- 
perature will cause quasilong-range order with a power 
law decay of correlation functions. If only one of the 
two amplitudes |$2;.y| has nonzero expectation value, 
the density Eq. (|2.1|) describes unidirectional (striped) 
CDW order. Again, the presence of a commensurate 
lattice potential makes the order long ranged at finite 
temperatures, whereas in the incommensurate case it is 
quasilong-ranged. 

In the incommensurate phase, in the absence of dis- 
order, the most general free energy density expanded in 
powers of and its gradients consistent with the sym- 
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FIG. 1: The homogeneous mean field solutions of Eq. (12.211 
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metries of the square lattice is given by 
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The homogeneous mean field solution of this model is 
summarized in Fig. ^ where the checkerboard, stripe and 
liquid phase values of ^x,y are shown along with the ac- 
companying free energy densities. 

In this work we focus on the influence of quenched dis- 
order on CDW order and thus consider adding a term to 
the free energy consisting of two complex random fields, 
coupling directly to ^x,y 



+ c.c. 



resulting in the total action 



(2.3) 



(2.4) 



where the complex random fields (/x = x, y) are pa- 
rameterized as H^(r) — /i^(r)e"'''(''\ are Gaussian 
distributed random variables with mean 0, and standard 
deviation and r;^ are uniformly distributed random 
phases on [0, 2n). 

Let us confine ourselves to the condensed phase where 
s < 0, and u > for stability. The interesting physics are 



encapsulated by the effects of altering the coupling con- 
stant V and the variance of the random field h^. While v 
changes the low energy ground states from checkerboard 
like configurations for w < {\^x\ = \^y\) to stripe like 
patterns (|$^| ^ \'^u\ — 0) for t; > 0, the variance 
should destabilize both types of states. 

A careful treatment of the coupling between CDW or- 
der and a random electrostatic potential yields random 
compression terms of the formi^ h^96^ (6^ = arg[<i>^]) 
omitted in the free energy Eq. (|2.3|l . In addition, an RG 
analysis of the full action JF^ + generates random 
shear terms. Random compression and shear terms are 
responsible for the power law decay of correlation func- 
tions on intermediate length scales, on which the influ- 
ence of topological excitations (dislocations) in the phase 
fields can be neglected^'*. 

In STM experiments, the correlation length of charge 
order is found to have values ranging iroxrX 2.5, to 
roughly'^ 5 CDW periods. In neutron scattering experi- 
ments, peak widths corresponding to correlation lengths 
larger then ten CDW periods were observed^. The cor- 
relation length describes the scale on which dislocations 
proliferate, and the presence of a relatively short correla- 
tion length indicates that there is no intermediate length 
scale on which compression and shear terms are impor- 
tant. For this reason, the omission of these terms from 
the elastic energy should be justified. 



III. NUMERICAL MINIMIZATION 

Due to the presence of the two complex random fields 
H^(r), we have elected to minimize Eq. H2.4|l numerically. 
The interplay of elastic and disorder energy causes frus- 
tration and gives rise to an exponentially large number of 
low lying states with similar energies but very different 
configurations. As these states are separated by large en- 
ergy barriers, relaxation after an external perturbation is 
very slow and glassy dynamics can be observed. For these 
reasons, numerically finding the ground state of such a 
system is a hard problem and as novel algorithms are de- 
veloped and employed, even to relatively simple models, 
new states with lower energies are inevitably foundiSi. 

Gradient methods which move strictly downhill in the 
energy landscape are fast, but prone to becoming stuck 
in local minima and are not always able to reproduce the 
results of slower ergodic methods. Simulated annealing 
algorithmsii have been the most successful at thoroughly 
sampling the possible configuration space by using a fic- 
titious temperature. By successively lowering this tem- 
perature, the resolution of finer and finer energy scales 
becomes possible while avoiding the danger of being stuck 
in a metastable excited state. 

As a compromise we have elected to employ a combi- 
nation of both greedy conjugate gradient and ergodic 
simulated annealing^^ methods. We allow for the possi- 
bility of local uphill moves where the configuration up- 
date involves making a downhill step in a random area 
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of the sample. The size of the randomly chosen region 
is annealed by tracking metropolis acceptance rate. This 
method faithfully reproduces quite quickly, the results of 
early Monte Carlo work on the random field XY model^". 

We have performed simulations for a number of lat- 
tice sizes L = {20,32,48,64,100}, commensurate with 
the experimentally observedi*Sii4 period of modulations 
in the local density of states of four lattice spacings 
{p = 4) and multiple realizations of disorder Nrd{L) = 
{200,200,150,150,100}. 

Let us consider aLxL square lattice of N sites labelled 
by i, then after rescaling to give dimensionless coupling 
parameters the continuum action Eq. (|2.4() in units where 
the lattice constant ao is set to unity and C3 — takes 
the form (with 9p,i = arg[$^,i]) 



fi—x^y I 



- 2/i^,i|$^,i| cos(e^^i - rif,^,) 



(3.1) 



with {i,j) indicating the usual sum over nearest neigh- 
bors and the factor of 1/2 is inserted to avoid double 
counting. The coupling matrix J^j- has diagonal el- 
ements Jx,i±x = Jy,i±y — Ci ^-nd off diagonal cou- 
plings Jx.i±y = Jy,i±x — C2- Wc havc chosen to set 
Ci = C2 = 1 and s = —0.1 thus restoring full rotational 
symmetry of the elastic energy on scales much larger than 
the lattice spacing and confining our analysis to the con- 
densed phase. We have also elected to pick the value of 
the quartic coupling u to ensure the condensation energy 
remains constant across the critical line u = by setting 
u{v > 0) = -s and u{v < 0) = -(s + v/2). 



IV. RESULTS 

Employing this minimization procedure we obtain sta- 
ble low energy field configurations like the ones shown 
in Figs. latoH For v = -0.1 and /iq = 0.6 (Fig. 12) 
we observe small circular regions less then 5 lattice spac- 
ings across, where the local random field configuration 
has taken a value that makes it energetically favorable to 
suppress the magnitude of either <^x or $y in its vicin- 
ity. For the uncoupled theory (Fig. ISJ '^x,y appear much 
smoother and although there still exists separate regions 
of stripe and checkerboard order, the interfaces between 
these regions are poorly defined. We may also observe the 
effects of positive w = 0.1 at the same disorder strength 
as seen in Fig.0] In this case, the magnitude of either $^ 
or is suppressed to zero over large regions of the sam- 
ple leading to a ground state field configuration with well 
defined domains having unidirectional order in either the 
X- or y-directions. 
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FIG. 2: The ground state field configuration on a 100 x 100 
lattice for a random field standard deviation of ^0 = 0.6 and 
a;j/- interaction v = —0.1 for a particular realization of disor- 
der. 




FIG. 3: The ground state field configuration on a 100 x 100 
lattice for a random field standard deviation of ho = 0.6 and 
xy-interaction v = 0.0 for a particular realization of disorder. 



We can construct the form of the density fiuctuations 
5p{r) (Eq. (|2.1|) ) from the minimized value of the spa- 
tially dependent order parameters to determine the ef- 
fects of altering v and /iq as seen in Fig. [S] For u < 
and hf) = 0.0 we observe robust checkerboard ordering 
which is coherent over the entire lattice. As the variance 
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FIG. 4: The ground state field configuration on a 100 x 100 
lattice for a random field standard deviation of ho = 0.6 and 
a;i/-interaction v — 0.1 for a particular realization of disorder. 
The phase of ^x,y fluctuates strongly in regions where its 
amplitude is suppressed. 



of the random field is increased, dislocations in the phase 
of the fields gradually destroy the local ordering, and 
the correlation length is reduced to less than three CDW 
periods for ho = 1.2. For v — 0.0 and Hq — 0.0, the mean 
field critical value in the disorder free theory, density fiuc- 
tuations in the x- and y-direction have a period of four 
lattice spacings and identical amplitudes. However, the 
response of the system to increasing ho is very different 
than in the v = —0.1 case. For any finite amount of dis- 
order, there is a significant enhancement in the size of 
unidirectional correlations as the sample breaks up into 
regions with the magnitude of either or $y greatly 
reduced. In the clean limit with w > we obtain purely 
unidirectional density fluctuations. As the strength of 
disorder is increased, large domains of 7r/2 relatively ori- 
entated appear and finally for large values of ho the sys- 
tem breaks up into many such regions of varying sizes. 

The presence of such strong stripeness leads to the nat- 
ural definition of a local Ising-like order parameter— 



E(r) 



|$.(r)|^-|$,(r)|^ 
|$,(r)|2 + |$j,(r)|2 



(4.1) 



which measures the tendency of the system to have only 
one of either or nonzero over some finite area, with 
a value between -1 and 1. 

As a result of our direct minimization procedure, we 
have calculated a full set of low energy field configura- 
tions for multiple system sizes and realizations of disor- 
der with zy-couplings v = {—0.1,0.0,0.1,0.2} and ran- 
dom field variances ho between 0.0 and 2.0. Using this 



information we can construct the disorder averaged cor- 
relation functions for each system size throughout the 
relevant parameter space. Two distinct types of correla- 
tions are of interest. The first are simple CDW correla- 
tion functions between the <i>p fields given by 



G^(0 - (<f>^(r)$;i(0)) 



(4.2) 



where fi € {x, y} while the second type measures fluctu- 
ations of the Ising-like order parameter Eq. H4.1|l 



Gs(r) = (E(r)S](0)) 



(4.3) 



with (■ • ■ ) indicating a spatial average in the ground state 
and the over-line denotes an average over multiple real- 
izations of disorder Nrd{L). As L becomes large both 
($^(r)} and (S(r)) tend to zero, and the connected and 
disconnected correlation functions are equivalent. Note 
that the definition of Gs (r) distinguishes between regions 
with strong unidirectional order with 7r/2 relative orien- 
tation. 

For sufficiently large variances in the magnitude of the 
random field, we expect that both Eqs. H4.2|l and H4.3|l 
will be characterized by exponential decays of the form 



Gs(r) 



(4.4) 



where and are their respective correlation lengths. 

(r) and Gs (r) are shown for distances up to 40 lattice 
spacings for v = —0.1,0.0 and 0.1 in Fig. El 

All correlations appear to decay exponentially for ho > 
0.4, and the most striking differences between G^{r) and 

Gs{r) can be seen at r = by comparing (<&2(0)) and 

(S2(0)). The xy-interaction parameter v has little ef- 
fect on the scale of the background CDW order, while 
the background Ising-like order, measured by Gs(0), in- 
creases by three orders of magnitude as v changes from 
—0.1 to 0.1. In addition it appears that (after proper 
finite size scaling) G^(0) is essentially a monotonically 
decreasing function of ho from the checkerboard to stripe 
parameter regime. The two insets in Fig. clearly show 
however, that the slope of Gs(0) vs. ho changes from 
positive to negative near v = 0.0. 

In order to determine the decay constants associated 
with G^(r) and Gs{r) we have performed a discrete 
Fourier transform of the disorder averaged correlation 
functions and fit them to a Lorentzian in fc-space for 
each linear system size, xy— interaction and random field 
variance. The actual correlation length is assumed to be 
equal to the width of the Lorentzian. Finite size scal- 
ing was then performed for each value of v and ho, as 
shown for v = —0.1 and ho = 0.6 in Fig. d to extract 
approximate infinite system size values of and • 

The results of the finite size scaling procedure can be 
seen in Fig. [S] where we plot the decay constants asso- 
ciated with Gx and Gs as a function of random field 
variance ho for various values of v, the error-bars in the 
fits are on the order of the symbol sizes. 
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FIG. 5: The density fluctuations Sp{r) on a 100 x 100 lattice for v = —0.1 (left column), v — 0.0 (center column) and 
V = 0.1 (right column) with ho = 0.0, 0.6 and 1.2 from top to bottom for one random field configuration. The central row was 
constructed using the values of ^x,y from Figs. |5|through2] 



For a fixed value of v, both and decrease mono- 
tonically as a function of disorder strength. The de- 
pendence of correlation lengths on the xy couphng v 
for a fixed random field strength is more interesting, as 
both correlation lengths are non-monotonic functions of 
V. Changing v from —0.1 to 0.0, the correlation length 
increases by almost ten lattice spacings in the regime of 



moderate disorder strength. This increase can probably 
be attributed to the fact that for v — 0.0, our model de- 
composes into two completely uncoupled unidirectional 
CDWs with ordering in the x and y direction, respec- 
tively. As disorder has a weaker influence on a unidirec- 
tional CDW as compared to a checkerboard CDW, this 
dependence of the correlation length on the value of v is 
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FIG. 6: Selected and Ge correlation functions for v = —0.1,0.0,0.1 and ho = 0.4,0.6,0.8,1.0,1.2 (legend applies to all 
panels) which were averaged over 100 realizations of disorder in a 100 x 100 sample. The insets show (E2(0)) vs. ho for v = —0.1 
(upper right) and v — 0.1 (lower right). 



plausible. A further increase of v to positive values leads 
to a reduction of . 

As V controls the degree of stripe order, it is expected 
that the correlation length for the Ising-like order pa- 
rameter increases from only a few lattice spacings for 



V = —0.1 to 20 to 40 lattice spacings for v = 0.1 in the 
regime of moderate disorder strength. It is important to 
note that for v — —0.1 and v = 0.0, Ising correlations 
are exclusively due to disorder fluctuations. When v ap- 
proaches zero coming from negative values, fluctuations 



7 



100 




0.04 0.06 

1 /L 



0.10 



FIG. 7: Finite size scaling of the correlation lengths S^x, ^nd 
using system sizes of L = {20, 32, 48, 64, 100} for u = -0.1 
and ho — 0.6. y-intercepts were extracted to approximate 
L oo. 
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at w = 0.1 (see inset of Fig. ISj). The fact that these 
peaks occur at different values of v leads to the interest- 
ing situation that > for positive v. Hence, with 
respect to the analysis of experimental data > S,x is 
a clear signature for a striped charge order if no random 
potential was present. In the limit of large v, the system 
breaks up into domains with either or $j, order and 
the correlation lengths £,x and become equal to each 
other. 



EMPIRICAL DETERMINATION OF 
STRIPENESS 



Due to the definition of our effective model for density 
fluctuations on the square lattice we have the direct abil- 
ity to measure the value of Eq. (|4.1|l using the ground 
state values of the independent CDW order parameters 
^x,y However, as we wish to make contact with the STM 
measurements discussed in Section it would be useful 
to determine a method of analyzing Sp(r) directly which 
might expose any underlying local Ising-like correlations 
that are not readily apparent either in real or Fourier 
space. 

With this goal in mind, we define an effective local 
Ising-like order parameter S(r) at each point in the sam- 
ple through the following procedure with more detail pro- 
vided in an appendix. 

1. Surround each lattice site r with a 4 x 4 box Dj. 
(depicted in Fig. Ill|l which is "centered" at the relative 
(1, 1) point. 

2. Define a local density-density correlation func- 
tion which lives in each 4x4 box (the smallest that 
contains enough information to resolve the wavevectors 
{±Kx,±Ky}) that is arbitrarily assigned to the (1,1) 
point. The result is 4 x 4 matrices 



FIG. 8: The correlation lengths associated with the exponen- 
tial decay otGx{r) and G^{r) extracted from finite size scaled 
Lorentzian widths in Fourier space. Color and symbol type 
denote different values of v, closed symbols with solid lines 
refer to and open symbols with dashed lines correspond to 
^E- The inset shows the two types of correlation lengths plot- 
ted for fixed ho as a function of v. Note that the maximum 
values occur at different values of the xy coupling. 



of the two fields and <i>j, are increasingly independent, 
and Ising correlations increase. When moving further 
into the stripe regime, we find ^e(w = 0.1) > Cs(w = 0.2) 
for all values of /iq. This can be attributed to the sharp- 
ening of domain walls between stripey regions with rel- 
ative orientation 7r/2 allowing them to better accommo- 
date the value of the local random field, increasing their 
overall length. 

While ^x reaches a maximum value at v = 0, peaks 



Su{r, 



') = {Sp{r')5p{r"))^,_ 



'"eor 



(5.1) 



where {■ ■ ■)r' -r"eD^ indicates an average over all sites 
r' £ Uj. and whose rows and columns are labeled by 
the X and y components of r' — r" employing periodic 
boundary conditions for the box. The specific form of 
one of the matrices contributing to this sum is given in 
Eq. IHJ. 

3. Perform a local discrete Fourier transform of 
Eq. 15.1|l using only those points in 



r' en. 



Sa{ry)e 



— ih-: 



(5.2) 



for k e {±KxT^Ky} at each of the L x L boxes. 4. Fi- 
nally, define an effective local Ising-like order parameter 
as the difference of local structure factor amplitudes at 
±Kx and iKy scaled by their sum in each box 
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FIG. 9: A comparison between the direct (E, left-panels) and 
effective (E, right-panels) Ising-like order parameters for ha — 
1.0 and V = {—0.1,0.0,0.1} for a particular realization of 
disorder in a 100 x 100 sample. 



We may now directly compare Eqs. H4.1(l and H5.3|l for 
different values of v at fixed ho as seen in Fig. |51 
The similarity between S(r) and S(r) is striking and 
shows both environs of checkerboard order (dark re- 
gions, |<i>a;(r)| = |$j,(t')|) and stripe order (light regions, 
|$,(r)|^|$,(r)|=0). 

The agreement between the left and right panels in 
Fig.Elcan be further quantized by defining an equal point 
correlator 



Q 



S(r)I](r) - I](r) S(r) 



where as is the standard deviation of E(r) 



(5.4) 
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FIG. 10; A comparison of the disorder averaged magnitudes 
of E and E for a 100 x 100 sample as a function of the variance 
of the random field. Finite size scaling appears to have little 

effect on these results. Symbols show the value of (|E(r)|) 
while lines refer to (|E(r)|). 



Using this definition we find values for Q of 0.62, 0.92 and 
0.71 for ho = 1.0 and v = -0.1,0.0 and 0.1. The smaller 
correlations for v = ±0.1 are due to the fact that O^. 
cannot be constructed symmetrically about the site r as 
we need to resolve the specific wavevectors {±1^2;, iKy} 
and thus small regions with rapid changes in ^x.y (sharp 
domain walls) will be poorly reproduced by the effective 
field S. 

A further comparison can be made by examining the 
disorder averaged values of the magnitudes of the direct 
and effective Ising-like order parameters. 



\n = (|S(r)|) 
|S| = (|S(r)|) 



(5.6) 



which are shown in Fig. ^| for various values oi v as a. 
function of disorder. Again we observe significant agree- 
ment between the direct and effective Ising-Like order 
parameters, now having averaged over 100 realizations 
of disorder. At ho = 0.0 we recover the results that in 
the disorder free theory should be identically zero for 
?; < and equal to one for v > 0. Increasing disor- 
der causes a smooth increase in unidirectional order for 

V < and a reduction for w > with the effective order 
parameter having a slightly larger dependence on ho. For 

V = 0.0, the magnitude of Ising-like order quickly jumps 
to a value near 0.5. 

The concurrence between S(r) and the effective object 
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E(r) inferred from the scaled difference in local structure 
factor peaks amplitudes at wavevectors iK^ and iKy 
is not surprising in the context that they are both cal- 
culated from the same underlying complex fields ^x.y m 
the condensed phase. However, the utility of Eq. (|5.3() 
becomes apparent when considering the plethora of ex- 
perimental STM spectra where only the LDOS is mea- 
sured and the underlying order is a topic of hot debate. 
The current analysis of such spectra involves perform- 
ing a discrete Fourier transform over the entire field of 
view. In real materials, disorder plays an important role, 
and short density-density correlation lengths are often 
observed. Therefore, such a global Fourier transform dis- 
cards a large amount of relevant local information which 
could in principle be used to probe any hidden electronic 
structure. 



VI. THE UNCOUPLED THEORY 

In the limit of large disorder a number of results can 
be explained for the uncoupled theory with v = 0. 

Upon examination of the various correlation lengths in 
Fig. |H1 it is apparent that ^a, is almost twice as large as 
for V = 0. This can be understood with the help of an 
approximate analytical argument: for u = 0, the lattice 
model Eq. 1)3. l|l reduces to two uncoupled unidirectional 
CDWs in the x- and y-directions. Concentrating on the 
numerator of the Ising order parameter Eq. H4.1|l for the 
moment, in the case v — the correlation function of S is 
then proportional to (|$2:(r)p|$a:(0)|^). If fluctuations 
of were described by a simple Gaussian theory, this 
would imply = 2^s, which indeed is approximately 
found in Fig. |H1 

A similar approach can be used to account for the sat- 
uration of and to 0.5 for u = as seen in Fig. 1101 
Again if we assume uncoupled fluctuations in ^^.y de- 
scribed by a normal distribution with mean and vari- 
ance cr|, we can directly calculate the expectation value 
of S by performing integrals in polar coordinates 

= i (6 1) 



VII. CONCLUSIONS 

This paper has characterized the correlations in a dis- 
ordered CDW state on the square lattice as a function of 
the parameter v in J-^ (which measures the degree of uni- 
directionality of the CDW order, v < corresponding to 
checkerboard states), and the random field strength Hq. 
We introduced a number of diagnostics to study the na- 
ture of the disordered state: the correlation lengths, 
of the CDW order, the correlation length of the Ising 
order associated with the uni-directionality, the on-site 
amplitudes G^(0), Gs(0), of these orders, and also dis- 



cussed in Section Ivl how closely related quantities could 
be measured directly in experiments. 

Our results are presented in detail in Sections IIVI 
and0 and here we highlight some notable features: 

(i) As is clear from the insets of Fig. El for v < 0, the 
strength of the Ising order increases with random-field 
strength, while the opposite is true for v > 0. 

(ii) The correlation length > only for v > 0, and 
this can serve as a diagnostic for the sign of v in an anal- 
ysis of the experiments. 

{in) We showed that the correlator Sn in Eq. (|5.1|) can 
serve as a very faithful diagnostic of the structure of S, 
and this should easily enable us to place experimental 
observations in the appropriate parameter space of the 
present theory. 

A full interpretation of the experiments should away 
a direct analysis of the experimental data along the 
lines proposed here. Nevertheless, a comparison of the 
qualitative structure of the figures presented here (es- 
pecially Fig [SJ with e.g. the STM observations of 
Ca2-2;Naa;Cu02Cl2 by Hanaguri et al. does indicate that 
this system has a value of v close to zero and likely pos- 
itive. 

For the future, our approach offers an avenue to quan- 
titatively correlate the STM experiments with neutron 
scattering. In particular, after determining the appropri- 
ate parameter regime of from an analysis of the STM 
data, the resulting Sp{r) configurations can be used as an 
input to determining the dynamic spin structure factor, 
as discussed in recent works^^'^^. 

While this paper was being completed, we learned of 
related results obtained by Robertson et alm^. 
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APPENDIX A: LOCAL ISING-LIKE ORDER 
PARAMETER 

This appendix provides more detail on the method 
used in the calculation of the effective local Ising-like or- 
der parameter S(r) described in Sec. IVl 
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Consider Fig. II II which shows one of A'' 4 x 4 boxes 
where the sites contained within the box are given the 
relative labels through 15 and it is "centered" at the 
(1,1) point here labelled 5. 



FIG. 11: One of 4 x 4 boxes used in the calculation of For the particular case shown in Fig. [TT] with r' = r[„ 

Sn (r, r' - r"). Here, has relative site labels to 15 and ^j^g matrix which contributes to the average Sn (r, r'-r") 
is centered at site 5 with r' currently located at site 10. ^ ^3 written out explicitly as 



Sa{r,r' -r") = •■• + 
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/Spu + Spe 6p5 + Spr + 5pi3 + dpi5 2{6p4 + Spi2) 6p5 + Spr + Spis + Spi^ 
25p2 2{5pi+5p3) A5po 2{5pi+5p:i) 

5pii + 5p(i (5p5 + (5p7 + (5pi3 + (5/915 2{5pi + 5pi2) 5p5 + 5p'j + 5pi3, + 5piz 

\ 5pio Spg + Spii 2(5/98 Spg + Spn 



r 



where we have employed the shorthand notation Spi 



the box using Eq. H5.2|l to give Sn{r, k). This expression 



5p{r'^). After performing the average over all r' — r" e niay then be evaluated at and i-K'j, using Eq. 1)5. 3|) 

i.e. calculating all 16 matrices at each site, Sn (r, r' - r") to give the effective Ising-like order parameter S(r). 
is Fourier transformed over the relative coordinates in 
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